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O ■ ABSTRACT 

o 



We introduce a Monte Carlo model of nonlinear diffusive shock acceleration al- 
lowing for the generation of large-amplitude magnetic turbulence, i.e., AB ^ Bq, 
where Bq is the ambient magnetic field. The model is the first to include strong 



c/3 ■ wave generation, efficient particle acceleration to relativistic energies in non- 

relativistic shocks, and thermal particle injection in an internally self-consistent 



manner. We find that the upstream magnetic field Bq can be amplified by large 
factors and show that this amplification depends strongly on the ambient Alfven 
Mach number. We also show that in the nonlinear model large increases in B do 
not necessarily translate into a large increase in the maximum particle momentum 
a particular shock can produce, a consequence of high momentum particles dif- 
fusing in the shock precursor where the large amplified field converges to the low 
ambient value. To deal with the field growth rate in the regime of strong fluctua- 
tions, we extend to strong turbulence a parameterization that is consistent with 
the resonant quasi-linear growth rate in the weak turbulence limit. We believe 
our parameterization spans the maximum and minimum range of the fluctuation 
growth and, within these limits, we show that the nonlinear shock structure, 
acceleration efficiency, and thermal particle injection rates depend strongly on 
the yet to be determined details of wave growth in strongly turbulent fields. The 
most direct application of our results will be to estimate magnetic fields amphfied 
by strong cosmic-ray modified shocks in supernova remnants. 
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Subject headings: Supernova Remnants, cosmic rays, shock acceleration. X-ray 
and radio emission, MHD turbulence 

1. INTRODUCTION 

Recent observations and modeling of several young supernova remnants (SNRs) suggest 
the presence of magnetic fields at the forward shock (i.e., the outer blast wave) well in 
excess of what is expected from simple compression of the ambient circumstellar field, -Bcsm- 
These large fields are inferred from (i) spectral curvature in radio emission (e.g., Reynolds 
& Ellison 1992; Berezhko et al. 1999), (ii) broad-band fits of synchrotron emission between 
radio and non-thermal X-rays (e.g., Berezhko et al. 2003; Volk, Berezhko & Ksenofontov 
2005) (see also Cowsik & Sarkar 1980), and (iii) sharp X-ray edges (e.g., Vink & Laming 
2003; Bamba et al. 2003; Volk, Berezhko & Ksenofontov 2005; Ellison & Cassam-Chenai 
2005). While these methods are all indirect, fields greater than 500 fiG are inferred in 
Cas A and values of at least several 100 fiG are estimated in Tycho, Kepler, SN1006, and 
G347.3-0.5. If i?csM ~ 3 — 10 fiG, amplification factors of 100 or more may be required to 
explain the fields immediately behind the forward shocks and this is likely the result of a 
nonlinear amplification process associated with the efficient acceleration of cosmic-ray ions 
via diffusive shock acceleration (DSA). The magnetic field strength is a critical parameter in 
DSA and also strongly influences the synchrotron emission from shock accelerated electrons. 
Since shocks are expected to accelerate particles in diverse astrophysical environments and 
synchrotron emission is often an important emission process (e.g., radio jets), quantifying 
the magnetic field amplification has become an important problem in particle astrophysics 
and has relevance beyond cosmic-ray production in SNRs. 

If cosmic-ray (CR) production is as efficient as expected in theories of nonlinear DSA (e.g., 
Blandford & Eichler 1987; Jones & Ellison 1991; Malkov & Drury 2001), the CR pressure 
gradient in the shock precursor can do work on the incoming plasma and, in principle, place 
a large amount of energy in magnetic turbulence. In fact, if the DSA process is to work 
at all, magnetic turbulence must be self-generated on all resonant scale lengths to provide 
the scattering necessary to drive the particle distribution to isotropy. If the turbulence 
remains weak, i.e., if AB/B ^ 1, the self-generation of turbulence by CR streaming in the 
shock precursor can be adequately described using quasi-linear theory (e.g., Skilling 1975; 
McKenzie & Volk 1982). For large amplitude turbulence, however, this analytic description 
becomes questionable and less rigorous approximations must be made. 

In principle, a complete description of nonlinear DSA, including magnetic field amplifi- 
cation, is possible with particle-in-cell (PIC) simulations. In reality, however, these simula- 
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tions are still far too computationally demanding to produce realistic models of astrophysical 
sources such as supernova remnants. Approximate methods describing nonlinear DSA must 
be used. Here, we have developed a method which incorporates a phenomenological model of 
particle diffusion and turbulence growth (similar to Bell & Lucek 2001; Amato & Blasi 2006) 
in a fully nonlinear Monte Carlo model of DSA (e.g., Jones & Ellison 1991; Ellison, Baring & 
Jones 1996; Berezhko & Ellison 1999). While not a first-principles description of the plasma 
physics, the computational efficiency of our approach gives it a significant advantage over 
PIC simulations and semi-analytic techniques based on the diffusion approximation in that 
we can determine the shock structure, the injection of thermal particles, the acceleration of 
these particles to relativistic energies, and the magnetic turbulence, all self-consistently in a 
fully nonlinear steady-state model. 

An important advantage of the Monte Carlo simulation is that, in effect, it solves a more 
basic set of equations governing the shock structure and DSA than do techniques based 
on diffusion-convection equations (e.g., Ellison & Eichler 1984). There is no assumption of 
isotropy for particle distributions and this allows an internally self-consistent treatment of 
thermal particle injection. Injection is still phenomenological and depends on the assump- 
tions made for the particle pitch-angle scattering, but these assumptions are applied equally 
to all particles. Since there is only one set of scattering assumptions, the Monte Carlo 
technique eliminates a free "injection parameter" which is present in all models based on 
the diffusion approximation to set the injection efficiency. As we show, the strong feedback 
between injection, shock structure, and magnetic field amplification makes this property of 
the Monte Carlo technique particularly important. 

Our preliminary results show that magnetic field amplification results in an increase in 
the maximum CR momentum, Pmax? a given shock system can produce compared to the 
case without amplification. For the case we consider, acceleration truncated by a finite-size 
shock, this increase is not, however, as large as the increase in B in the shocked region 
since fields on scales of the upstream diffusion length of the highest energy particles strongly 
influence Pmax- Furthermore, the fact that B ranges smoothly from the unshocked value far 
upstream, Bq, to the amplified value, B2, downstream from the shock will have consequences 
for electrons since radiation losses will be greatest in the downstream region where particles 
spend a large fraction of their time.^ Thus, the maximum momentum protons and electrons 
obtain may be determined by two very different field strengths. 



^We consistently use the subscript for far upstream values, the subscript 2 for downstream values, and 
the subscript 1 for values immediately upstream from the subshock. Thus, the overall compression ratio is 
Ttot = uo/u2 and the subshock compression ratio is rsub = ui/u2, where u is the bulk flow speed. 
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We also show that the amphfication factor, B2/B0, increases with Alfven Mach number, 
that is, for a given shock speed and ambient density, a small Bq will be amplified more than a 
larger Bq. The ability to dramatically increase very low ambient fields may make it possible 
for reverse shocks in young SNRs to accelerate particles and produce relativistic electrons 
and radio synchrotron emission (see Ellison, Decourchelle & Ballet 2005, for a discussion of 
DSA at reverse shocks is SNRs). The injection and acceleration of electrons in amplified 
fields will be considered in subsequent work, here, we consider only proton acceleration. 

In this paper, we demonstrate the feasibility of applying Monte Carlo techniques to 
DSA with field amplification but have made a number of approximations regarding wave 
generation and have neglected wave damping. An important advantage of the Monte Carlo 
method is that it can be generalized to include more realistic descriptions of non-resonant 
wave growth, linear and nonlinear damping, and the calculation of the momentum and space 
dependent diffusion coefficient from the turbulent energy density, and work along these lines 
is in progress. While a more accurate description of the plasma physics will influence all 
aspects of the acceleration process, we expect that Pmax will be most strongly determined 
by the plasma physics details and thus will remain a critical problem for understanding the 
origin of galactic cosmic rays for some time. 



2. MODEL 

2.1. Assumptions for Magnetic Turbulence Generation 

Consider a steady-state collisionless shock propagating along a uniform component Bq 
of a stochastic magnetic field B. We only consider parallel shocks where Bq is along x and 
the shock face is perpendicular to Bq. As the unshocked plasma approaches the shock and 
experiences the pressure gradient in the cosmic-ray precursor, the energetic particles back- 
streaming from the shock cause fluctuations of the field, AB, to grow. In the linear regime, 
AB is perpendicular to Bq and the local rate of growth of energy in waves is proportional 
to the particle pressure gradient. The plasma motion associated with these field fluctuations 
is initially Alfvenic: transverse and incompressible, and it will remain as such as long as 

AB < Bq. 

As the perturbations grow and reach AB > Bq, however, it is likely that waves with wave 
vectors k not aligned with Bq will be generated, due to local CR pressure gradients along the 
total B = Bq + AB. With AB > -Bq, it becomes impossible to predict the average value of 
the transverse pressure gradients and the resulting magnetic field structure without knowing 
the relative phases of different wave harmonics. The problem is further complicated by the 
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fact that this longitudinal, compressible turbulence may produce a strong second-order Fermi 
particle acceleration effect which, in turn, can damp the longitudinal fluctuations (see, for 
example, Schlickeiser, Campeanu & Lerche 1993). 

These complications place a precise description of plasma turbulence beyond current 
analytic capabilities and lead to our goal of obtaining a realistic approximation that includes 
the essential nonlinear effects from efficient DSA but still allows particle acceleration over a 
large dynamic range. In order to accomplish this, we consider two limiting cases. The first 
assumes there is no longitudinal turbulence, in which case the wave growth rate is determined 
by the Alfven speed in the un- amplified field Bq. This gives a lower bound to the growth 
rate. The upper limit assumes that the turbulence is isotropic, in which case the growth rate 
is determined by the Alfven speed in the much larger amplified field i?cfr (defined below). 
The real situation should lie between these two cases and while we consider these limits, we 
do not explicitly include second-order Fermi acceleration in our calculations. 

To get an initial estimate of the shock structure, we use a Monte Carlo simulation of 
diffusive shock acceleration, assuming that the motion of particles can be described as pitch- 
angle scattering in the Bohm limit. The Monte Carlo simulation provides the distribution 
of fast particles at all positions relative to the viscous subshock and allows for a calculation 
of the particle pressure gradient, dPp{x,p)/dx, which drives the amplification process in 
the shock precursor.^ Then, we use dPp{x,p)/dx to calculate the total energy density in 
magnetic turbulence, ?7tot( 

function of position across the shock. 

This calculation uses a semi-analytical expression for wavelength-dependent growth of 
magnetic turbulent energy density U{x,k), which is similar to the well-known equation for 
Alfven wave growth (e.g., equation B.8 in McKenzie & Volk 1982). The complete expression 
is presented below as equations (11) and (12). The turbulence amplification term in them, 
which represents the streaming instability is: 



where p is the particle momentum, Utot{x) = U{x, k)dk, p{k) is the momentum of parti- 
cles resonant with wavevector /c, and much of the complicated plasma physics is contained 
in Vgi an unknown growth-rate coefficient with dimensions of speed. Once the wave energy 
density is determined (as described below), it is used to calculate an "effective" amplified 



^It's important to note that while we speak of "thermal" and "superthermal" particles, and may refer 
to superthermal particles as cosmic rays, the Monte Carlo simulation makes no distinction between thermal 
and superthermal particles. The same assumptions are applied to particles regardless of their energy. 
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magnetic field, Beg{x) = Airlltotix) , and this field is then used to calculate the particle 
diffusion coefficient, D{x,p), as a function of position and momentum. 

For the growth of Alfven waves in quasi-linear theory, Vg = Va, where Va = Bq/ ^^/A7Fp{x) 
is the Alfven speed calculated with the un-amplified field and p(x) is the matter density 
at position x. As mentioned above, taking Vq = Va ignores the transverse gradients of 
pressure and only the component of the pressure gradient along the initial magnetic field Bq 
is accounted for. This choice of Vq provides a lower limit on the amplification rate and was 
used in Amato & Blasi (2006). If, on the contrary, we define Vq using the amplified field, 
i.e., Vg = Bef[{x)/ ^y4:7Tp{x), it reflects the situation where the growth rate is determined by 
the maximum gradient of Pp{x,p) along the fluctuating field lines. This provides an upper 
limit on the wave growth rate and was used in Bell & Lucek (2001). 

A similar argument applies to the resonance condition. In order to obtain a solution, 
we assume a resonance relation between k and p even though there is unlikely to be such 
a simple relation in strong turbulence. Again, we can define two limiting cases where the 
resonant wavelength for particles of momentum p lies between Xrcs{p) = '^'^rg{Bo,p) and 
Kes{p) = '^'^i^g{Bcs,p), where Vg = pc/{eB) is the gyroradius. The real situation should 
lie between the two extremes for Vg and Arcs- For this preliminary work, we set Xrcs{p) = 
2nrg{Bo,p), but vary Vg between the two limits, i.e., we introduce a parameter, < /aif < 1, 
such that 

/alf} , (2) 

and Vg varies linearly between Va (for /aif = 0) and Bcs/ A/47rp(x) (for /aif = 1). 

Finally, we assume a Bohm model for diffusion. The mean free path of a particle with 
momentum p at position x is taken to be equal to the gyroradius of this particle in the 
amplified field, i.e., X{x,p) = rg{x,p) = pc/[qBcs{x)], and the diffusion coefficient is then 
D{x,p) = Xv/3, where v is the particle speed. Closure of the model is provided by using the 
newly calculated diffusion coefficient in the Monte Carlo simulation, calculating a new shock 
structure and pressure gradient profile dPp{x,p) /dx, and iterating until mass, momentum, 
and energy fluxes are conserved throughout the shock. 

2.2. Monte Carlo Simulation 

The basic Monte Carlo (MC) simulation used here is described in detail in a number of 
papers including Jones & Ellison (1991) and Ellison, Jones & Reynolds (1990); Ellison, Jones 
& Baring (1999). The model calculates the shock structure and nonlinear particle spectrum 
self-consistently, including thermal particle injection. As mentioned above, injection can be 
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treated self-consistently in this model because the Monte Carlo technique does not require 
the diffusion approximation and no distinction between thermal and superthermal particles 
is made. Until now, however, in parallel shocks we have assumed a spatially independent 
form for the diffusion coefficient, i.e., that the diffusion coefficient, D{p), is proportional to 
gyroradius: D{p) = Xv/3, where A = rjrgfi is the particle mean free path, Vg^ = pc/{qBo) 
is the gyroradius, and 77 > 1 is an arbitrary parameter determining the "strength" of scat- 
tering.'^ In fact, we assume A and then individual particles pitch-angle scatter with some 
maximum scattering angle which is set to give the assumed mean free path (see Ellison, 
Jones & Reynolds 1990; Ellison & Double 2004, for full details of the pitch-angle scattering 
process) . 

The Monte Carlo simulation obtains the shock structure or bulk flow speed u{x) by 
iteration, ensuring that mass, momentum, and energy ffuxes are conserved across the shock. 
When DSA is efficient, the shock is modified by the back pressure of accelerated particles 
and the ffow speed, u{x), becomes a "smooth" function of x, i.e., an upstream precursor with 
Pp ~ pu"^ forms. In addition to shock smoothing, the overall compression ratio of the shock, 
rtot, will increase when the acceleration is efficient (e.g., Eichler 1984; Jones & Ellison 1991) 
and rtot is obtained by iteration as well. Below, we show how our code has been generalized 
to include the modification of the diffusion coefficient by the buildup of magnetic turbulence. 



2.3. Magnetic Field Amplification 

To calculate the effect of the pressure gradient of particles on magnetic field ffuctuations, 
we start with an equation similar to equation (B.8) in McKenzie & Volk (1982): 

d d d dPcr - 

ot OX ox ox 

Here Uw = {ABY/Att is the energy density of the waves (assuming that the kinetic energy 
density of shear wave motion equals the magnetic energy density), = {3u/2 — Va)Uw is the 
wave energy ffux (the 3/2 represents the sum of the Poynting flux and the ffux associated 
with the transverse motion of plasma in Alfven waves), = 11^/2 is the magnetic pressure 
of waves acting on the plasma flow, and L represents wave energy losses (or gains) due 



•^Parallel shocks are those where the shock normal is parallel to the magnetic field direction and oblique 
shocks are ones where the field makes some angle to the normal. The MC simulation has been generalized 
for plane, oblique shocks, in which case B and D vary with x as the strength and angle the magnetic field 
makes with the shock normal vary (e.g., Ellison, Baring & Jones 1996; Ellison, Jones & Baring 1999; Ellison 
& Double 2004). 
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to processes other than compression and amphfication by the considered instabihty, and 
dPcr/dx is the pressure gradient of cosmic rays exciting the Alfven waves. This equation 
describes the growth of the energy in Alfven waves through their instabihty in the presence 
of CR streaming, and it assumes that all the waves are moving upstream with respect to the 
plasma at speed Va, so Va,x = —Va- 

Now, following Bell & Lucek (2001), we separate the turbulence into downstream- and 
upstream-moving structures, and define the energy density of these structures per wavenum- 
ber interval as f/+(x, /c) and f/_(x, fc), respectively, so that the total energy density of tur- 
bulence is 

/•oo 

Utot= / [U.ix,k) + U+{x,k)]dk . (4) 
Jo 

We also define the partial pressure of particles with momentum p per unit momentum interval 
as Pp(yX^ p) so that the total pressure in particles, including thermal ones, is 



p,tot 



Pp{x,p)dp 



(5) 



To derive equations for U±{x,k), we apply a steady-state version of (3) to waves with 
wavenumber k in the interval Ak. For the energy density of these waves in the first order of 
Ak we substitute 

U^ = [U+{x,k) + U-{x,k)]Ak . (6) 

For the energy flux, we write 

^u{x) - Vg^ U.{x, k) + (^-u{x) + U+{x, k) Ak , (7) 

and for the pressure of particles interacting with these waves, 

P,, = Pp{x,p)Ap . (8) 

Substituting these expressions into (3), ignoring the energy loss (gain) term L, then dividing 
both sides by Ak and taking the limit A A; ^ 0, we get 

1. (^u{x) - Vg^ U.ix, k) + (^-u{x) + Vc^ U+ix, k) 

M(a;)^ Qf^-(a;, k) + ^f/+(x, k)^ - v^^- 
where we have replaced Va,x with a weighted wave speed 

f/+(x, k) - f/_(x, k) 



d Pc,{x,p) 
dx 



dp 



dk 



(9) 



fwt(a;, k) = Vc 



G 



U+{x,k) + U_{x,k) 



(10) 
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in the driving term. Using logical extension of the definition of an "average" wave 

speed by Bell & Lucek (2001) applied to a narrow wavenumber range Ak. This, again, 
is justified as long as AB ^ Bq, but becomes less clear for strong turbulence. In these 
equations, Vq is defined as in equation (2). 

The coefficient \dp/dk\ that the driving term has acquired is necessary to relate the inter- 
val of wavenumbers Ak of amplified waves to the interval of momenta of particles interacting 
with these waves Ap. As mentioned above, the relationship between p and k (i.e., the res- 
onant condition) is assumed to he k = l/r^o, where Vg^ = cp/^cBq), and Bq is the far 
upstream magnetic field. One may argue that the field that the particles 'feel' as uniform 
in the non-linear regime AB ^ Bq is the field carried by all long-wavelength (relative to 
the particle gyroradius) harmonics, rather than Bq. While this may be true in general, our 
present model is insensitive to the choice of the resonance condition due to the simplified 
"Bohm" form of the particle mean free path we assume. 

Equation (9) expresses energy conservation and accounts for the amplification of tur- 
bulence by the streaming instability with the growth rate determined from kinetic theory. 
Furthermore, we assume, as did Bell & Lucek (2001), that interactions between the forward 
and backward moving waves drive them to isotropy on a time scale ~ rg^o/Vc- In order to 
account for this interaction, we write equation (9) as the sum of the following two equations 
for U±{x, k): 



\U[X] 



dx\2 ^ ' 



Vr 



G 



U. 



U+ + U. 



-Vr 



G- 



dPcr{x,p) 



dx 
uix) + Vg 



dp 



dP 

-Vr " 



'G- 



[x,p) 



dk 



dp 



(f/_ - [/-,) 



dk 



V 



+ — (f/_ - f/+) 



(12) 



U+ + U.'^ dx 

which are solved iteratively in the MC simulation. We note that equations (11) and (12) are 
consistent with equation (9) with or without the relaxation terms on the right-hand sides of 
both equations, but these terms may become important in cases with small shock velocities. 

Equations (11) and (12) are generalizations of those introduced by Bell & Lucek (2001). 
The generalization has two essential improvements. First, it accounts for the spatial depen- 
dence of fiow speed u{x) due to nonlinear effects of efficient DSA and, consequently, treats 
'compression' of the amplified field adequately. To illustrate this effect, consider (11) and 
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(12), neglecting the cosmic-ray pressure gradient term and Va-, i.e., taking Vg ^ u. Adding 
the two equations then results in 

which can be easily integrated to give + oc n"'^/^, or i?cfr oc Consequently, in 

a shock with a total compression ratio rtot = 10, for example, the stochastic magnetic field 

3/4 

gets a boost in amplification by a factor of about r^^^ ~ 6 solely through the compression 
of the plasma. This compressional effect is especially important at the subshock and the 
change in magnetic turbulence energy density across the subshock will influence the subshock 
compression ratio r^uh- This in turn will have a strong effect on the injection efficiency in 
the Monte Carlo model. 

Second, we generalized the equations to describe the whole spectrum of turbulence 
U±{x,k) rather than a single waveband with Ak = k. We solve this system with a finite- 
difference method, integrating from far upstream (x — > — oo) to x. The quantities u{x) and 
Pp{x,p) are obtained from the Monte Carlo simulation, as described in Section 2.2. For 
simplicity in this initial presentation of our model, and because the shocks we are mainly 
concerned with have high Alfven Mach numbers {va ^ u), we have neglected Vq compared 
to u in the first and second terms in equations (11) and (12) in the numerical results we 
present here.^ We do, however, account for the wave or scattering center speed relative to 
the bulk fiow speed in determining the energy change a particle receives as it scatters in the 
converging fiow. In each interaction, we replace u{x) — > u{x) +v^t, and since Wwt is generally 
negative in the upstream region (i.e., f/_ > [/+), a finite Wwt dampens particle acceleration. 
Downstream from the shock we take v^t = 0. 

The initial condition for our equations is the far upstream magnetic turbulence spectrum 
U{x ^ oo, k). We take 

U.{x^-oo,k) = U+{x^-oo,k) = [ A:,< A;<A;^ 

^ ' ' I 0, k<k^ OT k> km , ^ ' 

where the limits km and k^ are chosen to encompass the range between the inverse gyroradii of 
thermal particles and the most energetic particles in the shock, respectively. For concreteness 
we take a = 1, but none of our results depend in any substantial way on a, km, or kc- Using 
our definition of the effective, amplified magnetic field, 

POO 

Bl^{x) =A7i [U4x, k) + U+{x, k)]dk , (15) 



"'ill all of the examples in this paper, Vq < 0.2u{x) in the shock precursor and Vq < 0.5u2 in the 
post-shock region. 
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the normalization constant A in equation (14) is determined by requiring that 



/ [U-{-oo,k) + U+{-oo,k)]dk = 



(16) 



2.4. Particle Scattering 



Once the amplified magnetic field is determined from equation (15), the scattering mean 
free path is set equal to the gyroradius of the particle in this field, i.e.. 



Equation (17) is essentially the Bohm limit and, as such, is a crude approximation. In a 
strong coUisionless shock, modified by efficient DSA, a significant fraction of the energy is 
contained in high-momentum particles. These particles have long mean free paths and will 
resonantly produce turbulence where long-wavelength harmonics contain most of the wave 
energy. Consequently, to low-momentum particles, the strong long- wavelength turbulence 
appears approximately as a uniform field and equation (17) is justified. For the highest 
energy particles, however, equation (17) will overestimate the scattering strength since, for 
these particles, most of the harmonics of the magnetic field appear as short-scale fiuctuations 
which are not very efficient at changing the particle's momentum. In this case, there is no 
reason to assume that the scattering is resonant, i.e., there is no simple resonance relation 
between k and p. This is a critical point since the form for X{x,p) at high p determines the 
maximum momentum that can be produced in a given shock system, and this is one of the 
important unsolved problems for DSA. 

Clearly, more physically realistic models for both the diffusion coefficient and the reso- 
nance condition are required for future work. An approach for determining the mean free 
path of particles in strongly turbulent fields is described in Bykov & Toptygin (1992), where 
non-resonant scattering and diffusive transport of particles in large-scale fiuctuations are 
taken into account. The model presented here, where the calculation of the power spectrum 
of turbulence U±{x,k) is coupled to the nonlinear shock structure, will be generalized to 
include more physically realistic wave-particle interactions in future work. 





(17) 



2.5. Momentum and Energy Conservation 



The total energy density of the MHD turbulence is defined by equation (4), and it 
is assumed to be equally shared between the stochastic magnetic field and the stochastic 
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incompressible transverse motion of the gas. The momentum flux of the turbulence is then 
the magnetic pressure, i.e., 



Pw,tot(x) = i ^ [U.{x, k) + U+{x, k)]dk = ^U, 



tot (3;) 



The total energy flux in turbulence is 

'3 



-^w,tot [x) 



-U 



or 



G 



[X) 



U4x,k) + 

3 

^ -u{x)Uu 



-u 



+ VG]U+{x,k) 



dk , 



[X) 



(19) 
(20) 



in the limit Vg{x) <^ u{x). 



In order to determine the bulk velocity profile, m(x), consistent with the back-reaction 
of accelerated particles and turbulence on the flow, the Monte Carlo simulation solves the 
following equations expressing the conservation of mass and momentum fluxes, i.e., 

p{x)u{x) = poUq, (21) 



p{x)u{xf + Pp,tot(a;) + 



Pouo , 



(22) 



where Pp,tot(a;) is the pressure produced by all particles, thermal and superthermal and Pq is 
the far upstream momentum flux. The particle pressure is related to Pp{x, p) by equation (5), 
and is calculated in the simulation directly from the trajectories of individual particles. The 
energy flux conservation relation is 

p{x)u{x)^ 



+ Fp^tot{x) + -Fw,tot(a;) + gc, 



PoMo 



wO 



(23) 



where -Fp^tot(a;) is the energy flux in all particles, gesc is the energy flux lost by particles 
leaving the system at the upstream free escape boundary (FEB), and Fq is the far upstream 
energy flux. In parallel shocks, u{x) can be determined from equations (21) and (22) alone. 
The relation between rtot and gesc (i-e., eq. 10 in Ellison, Moebius &: Paschmann 1990) allows 
equation (23) to be used to check the consistency of the simulation results, as we show with 
the examples below. 



3. RESULTS 

In all of the following examples we set the shock speed Uq = 5000 km s^^, the unshocked 
proton number density UpQ = lcm~^, and the unshocked proton temperature Tq = 10^ K. 
For simplicity, the electron temperature is set to zero and the electron contribution to the 
jump conditions is ignored. With these parameters, the sonic Mach number Mg ~ 43 and 
the Alfven Mach number M^if ~ 2300(l/iG/5o). 
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3.1. With and Without Magnetic Field Amphfication 

In order to obtain a solution which conserves momentum and energy fluxes, both the 
shock structure, i.e., u{x), and the overall compression ratio, rtot, must be obtained self- 
consistently. Figure 1 shows results for four shocks where u{x) and rtot have been determined 
with the iterative method described above. In all examples in this section, /aif = 0. First, we 
compare the results shown with heavy-weight solid curves to those shown with heavy-weight 
dotted curves. The heavy solid curves were determined with S-fleld ampliflcation while the 
dotted curves were determined with a constant Bes{x) = Bq. All other input parameters were 
the same for these two models, i.e., Bq = 30 fiG, /aif = (i.e., minimum wave ampliflcation), 
and an upstream free escape boundary at dpEB = — lO^rg(Mo), where rg(uo) = nipUoc/ (cBq). 
The most striking aspect of this comparison is the increase in Bes{x) when fleld ampliflcation 
is included (bottom panels). The magnetic fleld goes from Bes{x — oo) = 30 /iG, to 
i?efr > 1000 /xG for x > 0, and this factor of > 30 increase in B will influence the shock 
structure and the particle distributions in important ways. Note that there is about a factor 
of ~ 2 jump in B at the subshock at x = 0. As shown in Section 2.3, the jump in B at the 
subshock from compression (ignoring the contribution from the particle pressure gradient) 
is B2/B1 ~ {ui/u2)^^'^ = Tgy^, where ui is the flow speed immediately upstream from the 
subshock and the subshock compression ratio is deflned as rgub = Ui/u2- For the heavy- 
weight solid curve in Figure 1, rgub — 2.7 and B2/ Bi ~ 2, as observed. 

The solution without S-fleld ampliflcation (dotted curves) has a considerably larger r^t 
than the one with ampliflcation, i.e., for no i?-fleld ampliflcation, rtot — 22, and with 5-fleld 
ampliflcation (heavy solid curves), rtot — H.^ This difference in overall compression results 
because the wave pressure Pw.tot in equation (22) is much larger in the fleld amplifled case 
making the plasma less compressible. 

Two effects cause rtot to increase above the test-particle limit of 4 for strong shocks. 
The flrst is the production of relativistic particles which produce less pressure for a given 
energy density than non-relativistic particles making the plasma more compressible. The 
second, and most important, is the escape of energetic particles at the FEB. As indicated in 
the energy flux panels of Figure 1, the energy flux drops abruptly at x ~ c^feb as energetic 
particles diffuse past the FEB and escape the system. The energy lost from the escaping 
particles is analogous to radiation losses in radiative shocks and results in an increase in 
compression ratio. For high Mach number parallel shocks, 1/rtot ~ [5 — (9 + 16gcsc/-fo)^^^]/8 



^See Berezhko & Ellison (1999) for a discussion of how very large rtot's can result in high Mach number 
shocks if only adiabatic heating is included in the precursor. The uncertainty on the compression ratios for 
the examples in this paper is typically ±10%. 
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(EUison, Moebius & Paschmann 1990), so for the dotted curves rtot = 22 and gesc — O.SFq, 
and for the heavy sohd curves rtot = 11 and gcsc — O.S-Fq, consistent with the energy fluxes 
shown in Figure 1. 

The momentum flux is also conserved, but the escaping momentum flux is a much smaller 
fraction of the upstream value than for energy (Ellison 1985). Any departures from the far 
upstream momentum flux seen in Figure 1 are less than the statistical uncertainties in 
the simulation. As mentioned above, we only include adiabatic heating in the upstream 
precursor. If wave damping and heating were included, it is expected that rtot would be 
reduced from what we see with wave ampliflcation alone. 

For a given shock, an increase in rtot niust be accompanied by a decrease in the sub- 
shock compression ratio, rsub (see, for example, Berezhko & Ellison 1999, for a discussion 
of this effect). A large rtot nieans that high energy particles with long diffusion lengths get 
accelerated very efficiently and, therefore, the fraction of particles injected must decrease 
accordingly to conserve energy. The shock structure adjusts so weakened injection (i.e., a 
small rsub) just balances the more efficient acceleration produced by a large rtot. Since rsub 
largely determines the plasma heating, the more efficiently a shock accelerates particles caus- 
ing rtot to increase, the less efficiently the plasma is heated. For the cases shown in Figure 1, 
rsub — 3.8 for the amplified field case (heavy solid curves) and rsub — 2.4 for the case with 
no i?- field amplification (heavy dotted curves). 

In Figure 2 we show the phase space distributions, f{p), for the shocks shown in Figure 1. 
For the two cases with the same parameters except field amplification, we note that the 
amplified field case (heavy solid curve) obtains a higher Pmax and has a higher shocked 
temperature (indicated by the shift of the "thermal" peak and caused by the larger rsub) 
than the case with no field amplification (heavy dotted curves). It is significant that the 
increase in Pmax is modest even though B increases by more than a factor or 30 with field 
amplification. We emphasize that Pmax as such is not a parameter in this model; Pmax is 
determined self-consistently once the size of the shock system, i.e., cipEB, and the other 
environmental parameters are set. 

In order to show the effect of changing dpEB, "we include in Figs. 1 and 2 field am- 
plification shocks with the same parameters except that rfpEB is changed to — lOOOrg(Mo) 
(dashed curves) and —10^ rg{uo) (light-weight solid curves). From Figure 2, it's clear that 
Pmax scales approximately as dpEB and that the concave nature of f{p) is more pronounced 
for larger Pmax- The field amplification also increases with Pmax? but the increase between 
the dpEB = ~1000rg(Mo) and d^^B = — lO^rg(Mo) cases is less than a factor of two (bottom 
panels of Figure 1 ) . 
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In Figure 3 we show the energy density in magnetic turbulence, U+{x,k) + U-{x,k), 
the diffusion coefficient, D{x,p), and particle distributions as functions of k and p at three 
different positions in the shock. All of these plots are for the example shown with dashed 
curves in Figs. 1 and 2. The bottom panel shows how f{x,p) varies in the precursor where 
particles must diffuse upstream against the incoming plasma. The upstream diffusion length, 
[D{x,p)/u{x)]stvc, is some weighted average which is determined directly in the Monte Carlo 
simulation. The diffusion coefficients, shown in the middle panel, are determined from B^s 
(equation 15), where the U{x, /c)'s (top panel) are determined by solving equations (11) and 
(12) self-consistently with the shock structure using the particle pressure gradients deter- 
mined from f{x,p). The decrease in by ~ 40 between far upstream and downstream from 
the shock corresponds to the increase in Beg shown with the dashed curve in the bottom 
panel of Figure 1. The top panel clearly shows the spread in k where resonant interac- 
tions produce wave growth at the different x-positions, including the portion from thermal 
particles at high k. 

The efficiency of the shock acceleration process is shown in Figure 4. The curves on the 
left give the number density of particles with momentum greater than p, i.e., A^(> p), and the 
curves on the right give the energy density in particles with momentum greater than p, i.e., 
E{> p). Both sets of curves are determined solely from the downstream particle distributions 
(calculated in the shock reference frame) shown in Figs. 1 and 2 with heavy solid and 
dotted curves. Thus they ignore escaping particles and the energy in magnetic turbulence. 
Nevertheless, these curves indicate that the shocks are extremely efficient accelerators with 
> 50% of the energy density in f{p) placed in relativistic particles (i.e., p > rripc). The 
actual energy efficiencies are considerably higher since the escaping particles carry away a 
larger fraction of the total energy than is placed in magnetic turbulence. With gesc included, 
well over 50% of the total shock energy is placed in relativistic particles. Despite this 
high energy efficiency, the fraction of total particles that become relativistic is small, i.e., 
N{> p = rripc) ~ 10~^ in both cases. 

The effect of magnetic field amplification on the number of particles injected is evident 
in the left-hand curves. The larger rsub (solid curve) results in more downstream particles 
being injected into the Fermi mechanism with amplification than without. While it is hard 
to see from Figure 4, when the escaping energy flux is included, the shock with i?-field 
amplification puts a considerably smaller fraction of energy in relativistic particles than the 
shock without amplification. Again, injection depends in a nonlinear fashion on the shock 
parameters and the subshock strength will adjust to ensure that just the right amount of 
injection occurs so that momentum and energy are conserved. 
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3.2. Alfven Mach Number Dependence 

In Figure 5 we show three examples where Bq has been varied from 0.3 to 3 to 30 fiG; 
all other input parameters are kept constant including the physical distance to the FEB 
and /aif = 0. For these examples, IcipEBl = 1-7 x 10^° m = 5.6 x 10"'' pc. In units of 
rg(uo) = rripUQ/ {eBo), the units used for the x-coordinates in Figs. 1 and 5, this corresponds 
to I^febI = 1000, 100, and lOrg(Mo), for Bq = 0.3, 3, and 30 /iG, respectively. The top four 
panels showing u{x)/uq and energy flux have the same format as the corresponding panels in 
Figure 1. As in Figure 1, gesc is significant and rtot > 7 in all cases. The magnetic field panels 
differ from Figure 1 in that here we plot B{x)/Bq and it's clear that the amplification of B 
is greatest for the lowest Bq, i.e., B{x)/Bq increases with increasing Maif. For the examples 
shown here, B2/B0 ~ 400 for Bq = 0.3 /xG, B2/B0 ~ 150 for Bq = 3 /xG, and B2/B0 ~ 30 
for Bq = 30 /iG. In the bottom panels, we show the pressure in magnetic turbulence, Pw.tot, 
divided by the total far upstream momentum flux, Pq. For these examples, Pw,tot/-Po ^0.1 
and the magnetic pressure stays well below equipartition with the gas pressure. 

In Figure 6a we show the distribution functions corresponding to the shocks shown in 
Figure 5. As expected, the shock with the highest Bq = 30 /iG yields the highest Pmax- This 
Pmax, however, is only about a factor of 5 greater than that for Bq = 0.3 /iG obtained in 
a shock system of the same physical size. In Figure 6b we show the distribution functions 
for three cases where we have kept the upstream FEB boundary at the same number of 
gyroradii, i.e., (ipEB = — lOOrg(no). In this case, the situation with Pmax is reversed with the 
Bq = 0.3 /iG shock obtaining the highest Pmax- This is a combination of the fact that the 
physical size of the shock system is largest and that the field amplification is greatest for 
Bq = 0.3 //G. The shock structure results are not shown but are similar to those in Figure 5. 

3.3. Wave Amplification factor, /aif 

All of the examples shown so far have used the minimum amplification factor /aif = 
(equation 2). We now investigate the effects of varying /aif between and 1 so that Vq varies 
between Va{x) and B^six)/ ^y4:7rp{x). The other shock parameters are the same as used for 
the dashed curves in Figure 1, i.e., uq = 5000km s~^, Bq = 30 /iG, and dpEB = — lOOOrg(Mo). 

Figure 7 shows u{x)/uq and Bc{{{x)/Bq for /aif = 0, 0.1, 0.5, and 1 as indicated. The 
top panels show that increasing the growth rate (increasing /aif and therefore Vq) produces 
a large change in the shock structure and causes the overall shock compression ratio, rtot, to 
decrease. The decrease in rtot signifies a decrease in the acceleration efficiency and a decrease 
in the fraction of energy that escapes at the FEB, and the subshock compression adjusts 
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to ensure conservation of momentum and energy. The values for rtot and rsub are given in 
Table 1 and it is interesting to note that rsub increases as rtot decreases and becomes greater 
than 4 for /^if > 0.5. In contrast to the strong modification of u{x), there is little difference 
in Be{i{x)/BQ (bottom panels of Figure 7) and little change in Pmax (Figure 8), between these 
examples. 

The fact that increasing the wave growth rate decreases the acceleration efficiency shows 
the nonlinear nature of the wave generation process and comes about for two reasons. The 
most important is that the magnetic pressure term in equation (22), Pw,tot, becomes sig- 
nificant compared to p{x)u^{x) when /aif — > 1. The wave pressure causes the shock to be 
less compressive overall and forces rtot down. Any change in rtot changes the acceleration 
efficiency and therefore changes dPcr{x,p)/dx, the wave growth, and the shock structure. 
The second reason is that as the turbulence grows, w^t grows and the effective scattering 
center speed u{x) + fwt drops, causing the acceleration efficiency to drop. For the examples 
shown in Figure 7, the change produced by t>„t is modest compared to the effect from Pw.tot- 
Taken together, at least for the parameters used here, the increase in Pw,tot and fwt outweigh 
any increase in efficiency the amplified field produces. 

In Figure 8 we show the distribution functions for the four examples of Figure 7 and 
note that the low momentum peaks shift upward significantly with increasing /aif. The 
approximate momenta, measured in the downstream plasma frame, where the distributions 
peak, P / {mpCjp^^^, are listed in Table 1. As we have emphasized, the injection efficiency, i.e., 
the fraction of particles that enter the Fermi process, must adjust to conserve momentum 
and energy and the low momentum peaks shift as a result of this. The solid dots in Figure 8 
roughly indicate the injection point separating "thermal" and superthermal particles for the 
two extreme cases of /aif = and 1. The first thing to note is that this injection point is 
not well defined, a consequence of the fact that the MC model doesn't distinguish between 
"thermal" and "nonthermal" particles. Once the shock has become smooth, the injection 
process is smooth and the superthermal population smoothly emerges from the quasi-thermal 

Table 1: Effect of varying /aif. The errors on rtot and rgub are typically ±10%. 



/alf 


^tot 


^sub 


B2/B0 







9 


3.5 


30 


3.8 X 10^3 


0.1 


8 


3.7 


40 


4.4 X 10"3 


0.5 


6 


4.1 


50 


6.3 X 10-3 


1 


5 


4.3 


40 


8.0 X 10-3 
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population.^ Nevertheless, the approximate momentum where the superthermal population 
develops, Pinj, can be estimated and we mark this position with solid dots for /aif = 
and 1. What is illustrated by this is that the injection point shifts, relative to the post- 
shock distribution, when /^if is varied. This implies that, if injection is parameterized, the 
parameterization must somehow be connected to modifications in the shock structure. 

4. SUMMARY AND COMPARISON WITH ALTERNATIVE MODELS 

We have generalized a steady-state Monte Carlo model of efficient diffusive shock ac- 
celeration to include magnetic wave growth, allowing the wave amplitude to become large 
compared to the ambient field, i.e., AB/Bq ^ 1. The model uses a phenomenological treat- 
ment of wave generation by applying the linear growth rate formalism in the non-linear 
regime, but couples the nonlinear shock structure, the injection rate of thermal particles, 
the magnetic field amplification, and the determination of the maximum momentum par- 
ticles obtain from a physical constraint, i.e., the size of the shock system, in an internally 
self-consistent manner. 

Important limitations of our model are the simplified treatment of wave growth and the 
neglect of wave damping. We also assume a Bohm-like expression for the scattering mean 
free path [i.e., \{x,p) = pc/{eBcff)] rather than calculating D{x,p) from more fundamen- 
tal relations which account for particle motions in strongly anisotropic turbulence. More 
physically realistic models for turbulence generation and damping and particle scattering 
are certainly necessary and may well modify the results we present here. However, it is not 
straightforward to include these processes in nonlinear models and, to our knowledge, none 
have yet been presented with as extensive nonlinear coupling as we calculate. We believe the 
approximations we make are an important intermediate step and that the results we present 
are indicative of what more complete models will show. The fact that the MC technique can 
handle anisotropic particle distributions will be essential for including more precise plasma 
physics in future generalizations. 



^We note that the smooth emergence of a superthermal tail has been seen is spacecraft observations of the 
quasi-parallel Earth bow shock (i.e., Ellison, Mocbius & Paschmann 1990) and at interplanetary traveling 
shocks (i.e., Baring et al. 1997). 
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4.1. Particle-In-Cell Simulations 



The difference between our Monte Carlo simulation and a particle-in-cell (PIC) plasma 
simulation is that the pitch-angle diffusion is treated phenomenologically in the MC with an 
assumed scattering mean free path, rather than calculating particle trajectories in the self- 
generated magnetic field obtained directly from Maxwell's equations. This approximation is 
essential in order to make the MC technique fast enough to model acceleration over a wide 
dynamic range. We emphasize that while PIC simulations, in principle, can solve the shock 
acceleration problem completely, their computational requirements are extreme and they are 
not yet capable of modeling acceleration over a dynamic range large enough to model cosmic 
sources such as SNRs. The most important constraint on PIC simulations is that they must 
be done fully in three dimensions. As shown by Jokipii, Kota & Giacalone (1993) and Jones, 
Jokipii & Baring (1998), PIC simulations with one or more ignorable dimensions artificially 
confine particles to field lines and eliminate cross-field diffusion, an essential ingredient in 
diffusive shock acceleration, particularly in oblique shocks. The 3-D requirement means that 
all three box dimensions must be increased to accommodate high-energy particles with long 
diffusion length scales so that computing requirements become insurmountable for a large 
enough dynamic range. To model SNR shocks, they must be able to accelerate particles 
from thermal to highly relativistic energies and, if electrons are to be modeled, they must 
simultaneously include electron and proton scales.'' While PIC simulations will be able to 
investigate important problems, particularly those concerning injection, they will not be 
able to model the shock acceleration of electrons and protons to the energies necessary to 
produce broad-band radiation with parameters similar to those of SNRs in the foreseeable 
future. Until then, progress can be made with approximate methods. 



4.2. Semi- Analytic Models with B-field Amplification 

Besides the limited results from PIC simulations, the only models of magnetic field gen- 
eration with AB ^ Bo in nonlinear DSA that we are aware of are the semi-analytic results 
of Bell & Lucek (2001), Bell (2004, 2005), Ptuskin & Zirakashvih (2003), and Amato & Blasi 



^We note that the computational requirements for the relativistic shocks expected in 7-ray-bursts may 
actually be less stringent than those for the non-relativistic shocks in SNRs. In shocks with large Lorentz 
factors, particles start off relativistic (in the shock frame) and can gain a great deal of energy in just a 
few shock interactions. It is also possible that electron-positron plasmas dominate the 7-ray-burst fireball 
so that important results can be obtained without simultaneously covering electron-proton scales. In the 
non-relativistic shocks present in SNRs, however, both electrons and protons must be accelerated over a wide 
dynamic range from eV to TeV energies by crossing the shock many times. 
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(2006). As we have mentioned, our work uses the same basic wave generation formahsm as 
Bell & Lucek (2001), with an extensive generalization to include injection, nonlinear shock 
structure and particle distributions, and Pmax- 

More recently. Bell (2004, 2005) has attempted to improve the plasma physics by calcu- 
lating amplified fields from a so-called "directly driven" mode of wave instability.* In this 
scheme, upstream energetic particles standing in the shock frame produce a macroscopic 
current, and magnetic field, in the upstream plasma frame. This extended the Bell & Lucek 
(2001) analysis to wave modes other than resonant Alfvenic ones and emphasized that tur- 
bulence is likely the result of strongly driven, non-resonant modes at shorter wavelengths 
rather than Alfven waves. A recent application of this idea to relativistic shocks expected 
in GRBs is given Milosavljevic & Nakar (2006). While this is a promising idea for a non- 
resonant mechanism for wave amplification in strong shocks, careful consideration of the 
plasma return current must be made (e.g., F.C Jones, private communication). 

Another non-resonant mechanism for long-wavelength magnetic field fiuctuation ampli- 
fication in the vicinity of shock waves was studied by Bykov & Toptygin (2005). They 
demonstrated that a small fraction of neutral atoms can reduce the plasma transverse con- 
ductivity and result in a cosmic-ray current instability. This instability was shown to produce 
long-wavelength magnetic fiuctuations. 

An important attempt at describing turbulence with AB / B ^ 1 was made by Ptuskin 
& Zirakashvili (2003, 2005). These authors assumed a Kolmogorov-type nonlinear cascade 
and damping of self-generated turbulence by ion-neutral collisions. Most importantly, they 
obtained estimates for Pmax in a time- dependent analytic calculation. Nonlinear particle 
acceleration was assumed, but the analytic method required a number of approximations 
including the spatial distribution of energetic particles in the shock precursor and the spec- 
tral form of the energetic particles. Nevertheless, this work showed that nonlinear particle 
acceleration, combined with 5-field amplification, may strongly infiuence Pmax, one of the 
most important parameters in diffusive shock acceleration. 

In our estimation, the most highly developed model of nonlinear DSA with magnetic 
field amplification is that of Amato & Blasi (2006); work based on a series of papers by Blasi 
and co-workers (i.e., Blasi 2002, 2004; Blasi et al. 2005; Amato & Blasi 2005). ^ We now 



^Lucek & Bell (2000) and Bell (2004, 2005) also performed PIC simulations coupled to a 3-D MHD 
model of the background plasma. These results clearly showed that seed B-fields can be amplified by 
orders of magnitude but they were limited in dynamic range and did not self-consistcntly model the particle 
acceleration process. 

^We note that Amato & Blasi (2005, 2006) use a different technique from the previous Blasi et al. work 
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give a detailed comparison between that model and ours. The underlying assumptions for 
DSA are the same in both models, i.e., particles are driven to isotropy by interactions with 
magnetic turbulence in the background plasma. Particles gain energy when they diffuse 
in the converging flow and neither model includes second-order Fermi acceleration. Both 
models are also for plane-parallel, steady-state shocks. 

The most important difference revolves around the method of solution. We use a Monte 
Carlo technique while Blasi et al. solve a transport or diffusion-convection equation. In 
the Blasi et al. work, macroscopic quantities (pressure, energy flux, etc.) are derived as 
moments of the particle distribution function f{x,p), which is assumed to be isotropic in 
the shock reference frame. In contrast, the Monte Carlo simulation traces the stochastic 
motion of individual particles as they pitch-angle scatter off the background turbulence and 
calculates f{x,p) and its moments directly from the particle trajectories without making 
any assumptions about the isotropy of f{x,p). In effect, the MC simulation solves the more 
fundamental Boltzmann equation (e.g., Ellison & Eichler 1984). Of course, the semi-analytic 
technique is much faster computationally than the Monte Carlo technique and this will be 
important for building models in many applications. 

Injection can be treated self-consistently in the MC simulation because, once the pitch- 
angle scattering assumptions are made, they are applied equally to all particles and the 
number and energy of injected particles is fully determined. No distinction is made between 
thermal and superthermal particles and the viscous subshock is assumed to be transparent 
so there is a nonzero probability for any downstream particle with f > ^2 to be injected. 
The diffusion approximation, on the other hand, with its assumption of isotropy forces 
additional assumptions if injection is to be modeled. Blasi et al. treat the subshock as 
having a finite thickness comparable to a thermal particle's gyroradius and the injection rate 
is parameterized such that only those particles get injected, which have a gyroradius large 
enough to span the subshock. While this adds an additional parameter independent of the 
diffusion properties, it has been shown that parameters can easily be determined so that the 
Monte Carlo and semi-analytic models give similar results (see Ellison, Blasi & Gabici 2005). 

Another important difference is that Amato & Blasi (2006) have included a phenomeno- 
logical description of turbulent S-field heating, similar in implementation to that used in 
Berezhko & Ellison (1999), in addition to adiabatic heating. We only include adiabatic 
heating. The amplified turbulence may be dissipated through collisional and/or collisionless 
mechanisms and these include: (i) linear and nonlinear Landau damping (e.g., Akhieser et 



that allows an exact solution for an arbitrary choice of both the spatial and momentum dependent diffusion 
coefficient. 
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al. 1975; Achterberg & Blandford 1986; Kulsrud 1978; Vainshtein, Bykov & Toptygin 1993; 
Zirakashvili 2000), (ii) particle trapping (e.g., Medvedev 1999), and (iii) ion- neutral wave 
damping (e.g., Drury, Duffy & Kirk 1996; Bykov & Toptygin 2005). It's important to note 
that even though the damping rates tend to be smaller than the wave growth rates for long- 
wavelength fluctuations in turbulent plasmas with a substantial nonthermal component (e.g., 
Bykov & Toptygin 2005), the generation of magnetic energy density can be approximately 
balanced by the convection of that turbulence through the shock (as was assumed here) and 
MHD cascade processes. The heating of the precursor plasma by dissipation of small scale 
fluctuations modifies the subshock Mach number (e.g., Ellison, Berezhko & Baring 2000) 
and this in turn modifies injection. The overall acceleration efficiency and, of particular im- 
portance for X-ray observations, the temperature of the shocked plasma (e.g., Decourchelle, 
Ellison & Ballet 2000; Hughes, Rakowski & Decourchelle 2000; Ellison, Decourchelle & Bal- 
let 2004) will depend on wave dissipation. We plan to include physically realistic models of 
wave damping in the MC simulation in future work. 

The two most important elements in these nonlinear models is the turbulence generation 
and the diffusion coefficient that is derived from it. The assumptions for wave generation 
by the streaming instability are essentially identical in the Blasi et al. model and ours.^° 
The most important difference in the models lies in the calculation of the mean free path 
and diffusion coefficient. Here, we assume Bohm-like diffusion with Beg (equation 17), while 
Amato & Blasi (2006) use a mean free path due to resonant scattering that, in our notation, 
would be ^ 

where ^ ^ 

kres — — (25) 

rgfl cp 

is the resonant wavenumber. Both of these assumptions are approximations. Our Bohm-like 
expression is just a phenomenological recipe for particle diffusion in strong turbulence, while 
Ares is only formally valid for weak turbulence but is applied to diffusion in strong, anisotropic 
turbulence. An important step in advancing the state-of-the-art of nonlinear shock acceler- 
ation must center on improving the connection from wave generation to diffusion. 

That important aspects of DSA are sensitive to the assumptions made for diffusion can be 
seen by comparing the self-generated D{x, p) from these two models. Our use of equation (17) 



^°There is a minor difference in that Amato & Blasi (2006) assume that all waves generated by the 
streaming instability move upstream whereas we consider the interaction between waves moving upstream 
and downstream. 



- 23 - 



forces Bohm-like diffusion with D{x,p){x,p) oc vp and we show this dependence in Figure 3. 
Nevertheless, we have argued in Section 2.4 that equation (17) overestimates the scattering 
strength for the highest momentum particles because most of the harmonics appear as short- 
scale fluctuations to these particles. This would suggest that the actual p dependence of 
D{x,p) increases with p, likely as D{x,p) oc p"^ at the highest energies. In Amato & Blasi 
(2006), on the other hand, the momentum dependence of D{x,p) is shown to weaken at high 
p to the point where the scattering becomes strong enough that D{x, p) becomes independent 
of p in high sonic Mach number shocks. 

It is worth noting that the Amato & Blasi (2006) result is indicating that the magnetic 
turbulence spectrum has a slope approaching oc k~'^. For this turbulence spectral slope, 
the quasi-linear model predicts energy independent diffusion (in the weak fluctuation limit 
of course). Interestingly enough, the magnetic fluctuation spectrum we obtain with Bohm 
diffusion has, over a wide wavenumber range, a spectral slope roughly corresponding to oc 
(see Figure 3). So, if one formally calculated the quasi-linear diffusion from equation (24) 
with our fluctuation spectrum, the resulting diffusion coefficient would have only a weak 
energy dependence. This could mean that the magnetic fluctuation spectral shape is not 
very sensitive to the choice of diffusion model in the nonlinear calculations. Nevertheless, in 
this preliminary work, we prefer to use the Bohm diffusion model in the strong fluctuation 
limit. Since the behavior of D[x,p) at high p determines Pmax, as well as influencing all other 
aspects of the shock acceleration process, it's clear that much work remains to be done on 
this difficult problem. 

4.3. Other Models of Non-linear Shock Acceleration Without B-field 

Amplification 

There are a number of nonlinear models of DSA besides those discussed above. However, 
none of these models include -B-field amplification. 

In a series of papers, Berezhko and co-workers have applied a model of nonlinear DSA to 
broad-band observations of several young SNRs (e.g., Berezhko et al. 1996, 2002). They used 
a time- dependent solution of the cosmic-ray transport equation coupled to the gas dynamic 
equations in spherical symmetry and calculated the superthermal proton distribution, from 
the forward shock, at all positions in the remnant. From this they can determine the overall 
contribution a single SNR makes to the galactic cosmic-ray proton flux. Furthermore, since 
they included the acceleration of superthermal electrons, they calculated the photon emission 
from synchrotron and IC processes, as well as from proton-proton interactions and pion- 
decay. This allowed them to fit the broad-band photon observations and constrain the model 
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parameters. Of particular importance is their emphasis that magnetic fields much greater 
than typical ISM values are required to match broad-band photon observations. Equally 
important is their prediction, based on their fits, that TeV photon emission from several 
SNRs is most likely the result of pion-decay from protons rather than IC from electrons 
(e.g., Berezhko et al. 2003). 

This model has been extremely successful and has added to our understanding of SNRs, 
but it does make important approximations and simplifications. The model assumes Bohm 
diffusion, parameterizes the number of injected particles, and does not include wave ampli- 
fication.^^ 

Other nonlinear shock models include those of Kang and Jones and co-workers (e.g., 
Kang, Jones & Gieseler 2002; Kang & Jones 2006), Malkov and co-workers (e.g., Malkov 
1997; Malkov & Drury 2001; Malkov, Diamond & Jones 2002; Malkov & Diamond 2006), and 
the cosmic ray-hydrodynamical (CR-hydro) model of Ellison and co-workers (e.g., Ellison, 
Decourchelle & Ballet 2004). These models all involve different computational techniques, 
and all have their particular strengths and weaknesses. They have also been shown to 
produce similar nonlinear effects to those of Berezhko et al., and have been shown to be in 
quantitative agreement with our Monte Carlo model (before the addition of ^-amplification). 
Since none of these models yet include -B-field amplification we do not discuss them further. 



5. CONCLUSIONS 

We have introduced a model of diffusive shock acceleration which couples thermal particle 
injection, nonlinear shock structure, magnetic field amplification, and the self-consistent 
determination of the maximum particle momentum. This is a first step toward a more 
complete solution and, in this preliminary work, we make a number of approximations dealing 
mainly with the plasma physics of wave growth. Keeping in mind that our results are subject 
to the validity of our approximations, we reach a number of interesting conclusions. 



^^That i?-field amplification is not included in the Berezhko et al. model can easily be missed, particularly 
since the titles of some of these papers, e.g., "Confirmation of strong magnetic field amplification and nuclear 
cosmic ray acceleration in SN 1006" (Berezhko et al. 2003), might suggest that the model docs contain B-field 
amplification. The model of Berezhko et al. is a parallel shock model where the magnetic field is not explicitly 
included in the convection-diffusion equations. There is no i3-field amplification and the large downstream 
fields that Berezhko et al. infer from matching the observations are obtained by an ad hoc compression of the 
upstream field (see, for example, the discussion before equation (8) in Volk et al. 2002). To obtain 300 /iG 
downstream, for example, they must start with an unshocked field of ~ 50 fiG which is then compressed by 
the shock with a ratio typically rtot 6 for their models. 
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First, our calculations find that efficient shock acceleration can amplify ambient mag- 
netic fields by large factors and are generally consistent with the large fields believed to 
exist at blast waves in young SNRs, although we have not attempted a detailed fit to SNR 
observations in this paper. While the numerical values we obtain depend on the particular 
parameters for our examples, and we will investigate in detail how the amplification depends 
on sonic Mach number, age, and size of the shock system, etc., in future work, amplification 
factors of several 100 are clearly possible. 

More specifically, we find that the amplification, in terms of the downstream to far 
upstream field ratio B2/B0 is a strong function of Alfven Mach number, with weak ambient 
fields being amplified more than strong ones. For the range of examples shown in Figure 5, 
B2/B0 ~ 30 for Maif ~ 80 and B2/B0 ~ 400 for M^if ~ 8000. Qualitatively, a strong 
correlation between amplification and Maif should not depend strongly on our approximations 
and may have important consequences. In young SNRs, the expansion of the ejecta material 
will drastically reduce any original, pre-SN circumstellar magnetic field. In fact, for any 
conceivable progenitor, the magnetic field inside of the reverse shock will drop to values 
too low to support the acceleration of electrons to radio emitting energies only a few years 
after the explosion (e.g., Ellison, Decourchelle & Ballet 2005). Evidence for radio emission 
at reverse shocks in SNRs has been reported (see, Gotthelf et al. 2001, for example) and 
the strong amplification of low fields we see here, may make it possible for reverse shocks 
in young SNRs to accelerate electrons to relativistic energies and produce radio synchrotron 
emission. If similar effects occur in relativistic shocks, these large amplification factors will 
be critical for the internal shocks presumed to exist in 7-ray bursts (GRBs). Even if large 
5-field amplification is confined to non-relativistic shocks, which tend to be more efficient 
accelerators than relativistic ones, amplification will be important for understanding GRB 
afterglows since the expanding fireball will slow as it moves through the interstellar medium 
and will always go through trans-relativistic and non-relativistic phases. 

As expected, amplifying the magnetic field leads to a greater maximum particle momen- 
tum, Pmax, a given shock can produce. Quantifying Pmax is one of the outstanding problems 
in shock physics because of the difficulty in obtaining parameters for typical SNRs that allow 
the production of cosmic rays to energies at and above the CR knee near 10^^ eV. Assuming 
that acceleration is truncated by the size of the shock system, we determine Pmax from a 
physical constraint: the relevant parameter is the distance to the free escape boundary in 
diffusion lengths. This means that the limit on acceleration feeds back on the shock structure 
and also mimics, in terms of the spectral shape, what happens in actual shocks where f{x,p) 
must turn over smoothly at the highest energies (as in Figure 6). The spectral shape will be 
particularly important if the model is applied to the knee of the cosmic-ray spectrum or to 
nonthermal X-ray emission in SNRs. 
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Our results show that Pmax does increase when field amplification is included, but the 
increase is considerably less than the amplification factor at the shock B2/ Bq (compare 
the heavy dotted and heavy solid curves in Figure 2). The main reason for this is that 
high momentum particles have long diffusion lengths and the precursor magnetic field well 
upstream from the subshock strongly infiuences Pmax- We calculate the spatial structure of 
the amplified magnetic field (Figs. 1 and 5) and, as expected, field amplification is greatest 
near the subshock and -Bes merges into the ambient field far upstream. The diffusion length 
that determines Pmax (i-e., [-D(x,Pmax)/w(x)]avc) is some weighted average over the varying 
u{x) and Bcs{x) and is considerably greater than that estimated from B2 alone. If the shock 
size, in our case (ipEB; limits acceleration, Pmax will be considerably less than crude estimates 
using a spatially independent B2. 

On the other hand, particles spend a large fraction of their time downstream from the 
shock where the field is high and collision times are short. If shock age limits acceleration 
rather than size, we expect the increase in Pmax from the amplified field to be closer to the 
amplification factor, B2/B0. Thus, the spatial structure of the precursor field and u{x), in 
addition to the overall amplification of B, will determine the relative time spent upstream 
versus downstream and will determine Pmax for a given shock system. This will be even 
more critical for electrons than for ions since electrons experience synchrotron and inverse- 
Compton losses which will mainly occur downstream. Again, the qualitative nature of the 
above conclusions should not depend on the particular parameters and approximations we 
make here. We leave more detailed quantitative work for future study. 

Finally, it is well known that DSA is inherently efficient. Field amplification reduces 
the fraction of shock ram kinetic energy that is placed in relativistic particles but, at least 
for the limited examples we show here, the overall acceleration process remains extremely 
efficient. Even with large increases in Bef[{x), well over 50% of the shock energy can go into 
relativistic particles (Figure 4). As in all self-consistent calculations, the injection efficiency 
must adjust to conserve momentum and energy. In comparing shocks with and without field 
amplification, we find that field amplification lowers rtot and, therefore, individual energetic 
particles are, on average, accelerated less efficiently. In order to conserve momentum and 
energy, this means that more thermal particles must be injected when amplification occurs. 
The shock accomplishes this by establishing a strong subshock which not only injects a larger 
fraction of particles, but also more strongly heats the downstream plasma. This establishes 
a nonlinear connection between the field amplification, the production of cosmic rays, and 
the X-ray emission from the shocked heated plasma. 

The authors are grateful to P. Blasi and F.C Jones for useful discussions. DCE ac- 
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-27- 



(NNH04Zss001N-LTSA), and AB was supported, in part, by RBRF 06-02-16844. 

REFERENCES 

Achterberg, A. & Blandford, R. D. 1986, MNRAS, 218, 551 

Akhieser, A.I., Akhieser, I.A., Polovin, R.V., Sitenko, A.G. & Stepanov, K.N. 1975, "Plasma 
Electrodynamics", Pergamon Press, Oxford 

Amato, E. & Blasi, P. 2005, MNRAS, 364, L76 

Amato, E. & Blasi, P. 2006, MNRAS, in press, astro-ph/0606592 

Bamba, A., Yamazaki, R., Ueno, M. & Koyama, K. 2003, ApJ, 589, 827 

Baring, M. G., Ogilvie, K. W., Ellison, D. C. & Forsyth, R. J. 1997, ApJ, 476, 889 

Bell, A. R. 2004, MNRAS, 353, 550 

Bell, A. R. 2005, MNRAS, 358, 181 

Bell, A. R. & Lucek, S.G. 2001, MNRAS, 321, 433 

Berezhko, E.G. & Ellison, D.C. 1999, ApJ, 526, 385 

Berezhko, E. G., Elshin, V. K., & Ksenofontov, L. T. 1996, JETP, 82, 1 

Berezhko, E. G., Ksenofontov, L. T. & Petukhov, S. I. 1999, Proc. 26th Int. Cosmic- Ray 
Conf.(Salt Lake City), 4, 431 

Berezhko, E. G., Ksenofontov, L. T. & Volk, H. J. 2002, A&A, 395, 943 

Berezhko, E. G., Ksenofontov, L. T. & V51k, H. J. 2003, A&A, 412, 111 

Blandford, R. & Eichler, D. 1987, Physics Reports, 154, 1 

Blasi, P. 2002, Astroparticle Physics, 16, 429 

Blasi, P. 2004, Astroparticle Physics, 21, 45 

Blasi, P., Gabici, S. & Vannoni, G. 2005, MNRAS, 361, 907 

Bykov, A. M. & Toptygin, I. N. 1992, Sov.Phys. JETP, 74, 462 

Bykov, A. M. & Toptygin, I. N. 2005, Astronomy Letters, 31, 748 



- 28 - 

Cowsik, R. & Sarkar, S. 1980, MNRAS, 191, 855 
Decourchelle, A. Ellison, D.C. & Ballet, J. 2000, ApJ, 543, L57 
Drury, L. O'C, Duffy, R & Kirk, J. G. 1996, A&A, 309, 1002 
Eichler, D. 1984, ApJ, 277, 429 
Ellison, D.C. 1985, J. Geophys. Res., 90, 29 

Ellison, D. C., Baring, M. G. & Jones, F. C. 1996, ApJ, 473, 1029 
Ellison, D.C. Berezhko, E.G. & Baring, M.G. 2000, ApJ, 540, 292 

Ellison, D. C, Blasi, R & Gabici, S. 2005, proc. 29tli ICRC (Puna, India), astro-ph/0507107 

Ellison, D. C. & Cassam-Chenai, G. 2005, ApJ, 632, 920 

Ellison, D.C, Decourchelle, A. & Ballet, J. 2004, A&A, 413, 189 

Ellison, D.C, Decourchelle, A. & Ballet, J. 2005, A&A, 429, 569 

EUison, D. C. & Double, G. R. 2002, Astroparticle Rhysics, 18, 213 

Elhson, D. C. & Double, G. R 2004, Astroparticle Rhysics, 2, 323 

EUison, D. C. & Eichler, D. 1984, ApJ, 286, 691 

EUison, D. C, Jones, F. C. & Baring, M. G. 1999, ApJ, 512, 403 

EUison, D. C, Jones, F. C. & Reynolds, S. R 1990, ApJ, 360, 702 

EUison, D. C, Moebius, E. & Paschmann, G. 1990, ApJ, 352, 376 

Gotthelf, E. v., Koralesky, B., Rudnick, L., Jones, T. W., Hwang, U. & Petre, R. 2001, ApJ, 
552, L39 

Hughes, J. P., Rakowski, C. E. & Decourchelle, A. 2000, ApJ, 543, L61 
Jokipii, J. R., Kota, J. & Giacalone, J. 1993, Geophys. Res. Lett., 20, 1759 
Jones, F. C. & EUison, D. C. 1991, Space Science Reviews, 58, 259 
Jones, F. C, Jokipii, J. R. & Baring, M. G. 1998, ApJ, 509, 238 
Kang, H. & Jones, T. W. 2006, Astroparticle Physics, 25, 246 



- 29 - 



Kang, H., Jones, T. W. & Gieseler, U. D. J. 2002, ApJ, 579, 337 

Kulsrud, R. M. 1978, Astronomical Papers Dedicated to Bengt Stromgren, eds. Reiz, A. and 
Andersen, T. p. 317 

Lucek, S. G. & Bell, A. R. 2000, MNRAS, 314, 65 

Malkov, M. A. 1997, ApJ, 485, 638 

Malkov, M. A. & Diamond, R H. 2006, ApJ, 642, 244 

Malkov, M. A., Diamond, R H. & Jones, T. W. 2002, ApJ, 571, 856 

Malkov, M. A. & Drury, L.O'C. 2001, Reports of Progress in Physics, 64, 429 

McKenzie, J. F. & Volk, H. J. 1982, A&A, 116, 191 

Medvedev, M. V. 1999, Physics of Plasmas, 6, 2191 

Milosavljevic, M. & Nakar, E. 2006, submitted to ApJ, astro-ph/05 12548 
Ptuskin, V. S. & Zirakashvih, V. N. 2003, A&A, 403, 1 
Ptuskin, V. S. & Zirakashvih, V. N. 2005, A&A, 429, 755 
Reynolds, S. P. & Elhson, D. C. 1992, ApJ, 399, L75 
Schlickeiser, R., Campeanu, A. & Lerche, L. 1993, A&A, 276, 614 
SkiUing, J. 1975, MNRAS, 172, 557 

Volk, H. J., Berezhko, E. G., Ksenofontov, L. T. & Roweh, G. P. 2002, A&A, 409, 563 

Vainshtein, S. I., Bykov, A. M. & Toptygin, I. 1993, "Turbulence, current sheets, and shocks 
in cosmic plasma", Gordon and Breach Science Publishers, Langhorne, Pa., U.S.A. 

Vink, J. & Laming, J. M. 2003, ApJ, 584, 758 

Volk, H. J., Drury, L. O. & McKenzie, J. F. 1984, A&A, 130, 19 

V51k, H. J., Berezhko, E. G. & Ksenofontov, L. T. 2005, A&A, 433, 229 

Zirakashvih, V. N. 2000, JETP, 90, 810 



This preprint was prepared with the AAS I^TJt;X macros v5.2. 



- 30 - 



o 
3 



X 
13 



O 
CL 
\ 

X 

3 



E 
o 



0.5 




5 

2 
1 

0.5 

0.2 
0.1 

5 

2 
1 

Ld 

0.2 
0.1 

lO""' 

CO ' 
CO 

D 



X 

3 




hv. solid: B— amp, r^jj^=11 
dotted: no B— amp, rjQt=22 

— ^WW— Hi....... 



10" 



m 



10 



-5 



H \ ^ \ H J I I I I I I I 





I ' I ' I ' I 

subshock 



I I I I I I I I 







-6 -4 -2 



-4 -2 2 4 



Fig. 1. — Shock structure including momentum and energy fluxes (in units of far upstream 
values) and the effective magnetic field vs. x. Note that the horizontal scale has units of 
rg(Mo) = rripUo/ {cBq) and is divided at x = —5rg{uQ) between a linear and log scale. In all 
panels, the heavy dotted curves show results without amplification and all other curves are 
with amplification. The curves showing the energy flux drop sharply at the upstream FEB, 
which is at — lO'^rg(uo) for the heavy solid and dotted curves, at — lOOOrg(uo) for the dashed 
curves, and at — lO^rg(Mo) for the light solid curves, as particles freely leave the system. 
When this escaping flux is included, energy and momentum are conserved to within ±10%. 
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Fig. 2. — Phase space distributions for the shocks shown in Figure 1. These spectra are 
multiphed by [p/(?Tipc)]^ and are calculated downstream from the shock in the shock rest 
frame. As in Figure 1, the heavy solid and dotted curves have c^feb = — lO'^rg(Mo), the 
dashed curve has (ipEB = — lOOOrg(Mo), the light solid curve has dpEB = — lO^rg(Mo). 
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Fig. 3. — The top panel shows U+{x,k) + U-.{x,k) vs. k at three positions relative to 
the subshock. Here and in the other two panels, the solid curve is calculated downstream 
from the shock, the dashed curve is calculated at x = —rg^uo) upstream from the subshock, 
and the dotted curve is calculated at x = — lOOrg(uo) upstream from the subshock. The 
middle panel shows the diffusion coefficient with an additional dash-dotted curve showing 
the far upstream value. The bottom panel shows the distribution functions, multiplied by 
[p/{'mpc)]*, at the various positions. These distributions are calculated in the shock frame. 
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Fig. 4. — The left-hand curves show the fraction of particles with momentum greater than p, 
while the right-hand curves show the fraction of energy density in particles with momentum 
greater than p, for the shocks shown in Figs. 1 and 2 with heavy solid and dotted curves. 
These curves are calculated from the distributions shown in Figure 2 and do not include the 
particles that escaped at (ipEB- While the number fraction of escaping particles is small, the 
energy fraction is significant, as indicated in Figure 1. The solid curves show results with 
field amplification and the dotted curves show results without field amplification, both with 



-34- 





1 


o 








— 

X 













o 
Li_ 


5 


\ 

X 
D 


2 


Li_ 


1 


>^ 




Ol 
i_ 


0.5 


0) 




c 
LJ 


0.2 


1000 


o 






100 






X 






10 


CD 






1 


o 





CL 




\ 
o 


-2 


* 

Q_ 


-4 


o 






-6 


O 




_l 


-8 



. solid: Bq=0.3/^G 
. dashed: Bq=3//G 
, dotted: Bo=30//G 

. 1 I I I I I I I I I I I I I . 




I I I I I I I I I I I 



I I I I I I [ I I I I I r|_ I I [ I I I I I J 




J I I I I I I I I I I L. 



I I I I 




I.I 



— I 



-3 -2 -1 -4-2 2 
-Log^o (-x) X [r (uq)] 



Fig. 5. — Comparison of shocks with different far upstream fields Bq. In all panels, the solid 
curves are for Bq = 0.3 /iG, the dashed curves are for Bq = 3/iG, and the dotted curves 
are for Bq = 30 fiG. The FEB is placed at the same physical distance in all cases with 
c^FEB = —1.7 X 10^° m. The horizontal axis is in units of rg{uQ) = mpUo/ (cBq), and is split at 
X = — 5rg(Mo) between a linear and logarithmic scale. Note that Bes increases most strongly 
for Bq = 0.3 fiG, but that the pressure in magnetic turbulence never gets above ~ 10% of 
the total pressure. The overall compression ratios are: rtot — 9 for Bq = 0.3 /iG, rtot — 12 
for Bq = 3 fiG, rtot — 8 for Bq = 30 /iG, values consistent, within statistical errors, with gcsc, 
as indicated in the energy flux panels. 



- 35 - 




Fig. 6. — The top panel shows the distribution functions obtained in the shock frame in 
the downstream regions of the shocks shown in Figure 5. These shocks have the FEB at 
the same physical distance from the subshock and the shock with Bq = 30 fiG produces the 
greatest Pmax- In the bottom panel we show a similar set of curves only here the FEB was set 
at a fixed dpEB = — lOOrg(Mo) upstream. In this case, the shock with Bq = 0.3 //G produces 
the greatest Pmax, a result of the larger shock size and greater amplification factor the high 
Maif shock receives. 



- 36 - 




Fig. 7. — Shocks with varying f^u as indicated. In all cases, Uq = 5000km s ^, i?o = 30 /iG, 
and dpEB = — lOOOrg(uo). 
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Fig. 8. — Distribution functions for the shocks shown in Figure 7 calculated downstream 
from the shock in the shock reference frame. The solid dots give the approximate position 
for the transition between "thermal" and superthermal particles for the two extreme cases 
of /aif = and 1. 



